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Abstract. Scmidefinite programs are an important class of convex optimization prob- 
lems. It can be solved cfBciently by SDP solvers in Matlab, such as SeDuMi, SDPT3, 
DSDP. However, since we are running fixed precision SDP solvers in Matlab, for some 
applications, due to the numerical error, we can not get good results. SDPTools is a 
■ Maple package to solve SDP in high precision. We apply SDPTools to the certification 

^ of the global optimum of rational functions. For the Rumps Model Problem, we obtain 

the best numerical results so far. 



u 
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1. Preliminaries 

SDPTools is a Maple package to solve a semidefinite program (SDP): 

T 

c a 

s.t. F{x) > 0, 

where 

m 

F{x) = Fo+Y,^^P^■ 



> 

(N 
00 

. minimize x 



i=l 

nxn 



^ , The problem data are the vector c E M™ and m + 1 symmetric matrices Fq, . . . , Fm G 

H ' The dual problem associated with the semidefinite program (1) is 

maximize —TiFgZ 

s.t. TrFiZ = Ci, i = l,...,m, (2) 
Z >0. 

Here the variable is the matrix Z = Z^ E M"^", which is subject to m equality constraints 
and the matrix nonnegativity condition. 

The SDP (1) and its dual problem (2) can be solved efficiently by algorithms in Se- 
DuMi [9], SDPT3 [10], DSDP [2], and SDPA [3]. However, since we are running fixed 
precision SDP solvers in Matlab, we can only obtain numerical positive semidefinite ma- 
trices which satisfy equality or inequality constrains approximately. For some applications, 
such as Rump's model problem [7], due to the numerical error, the computed lower bounds 
can even be significantly larger than upper bounds, see, e.g.. Table 1 in [4]. These numerical 
problems motivate us to consider how to use symbolic computation tools such as Maple to 
obtain SDP solutions with high accuracy. 
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2. An algorithm for solving SDPs 

Semidefinite programs are an important class of convex optimization problem for which 
readily computable self-concordant barrier functions are known, so interior-point methods 
are applicable. Our algorithm for solving SDPs is based on the potential reduction methods 
mentioned in [11], which is to minimize the potential function (55) in [11] 

<p{x, Z) = vVnlog(TrF(x)Z) + ^{x, Z) (3) 
= (n + v\/n)\ogiTrF{x)Z) — log detF(x) — log detZ — nlogn 

The first term in the right hand of the first equation is the duality gap of the pair x, Z 
and 'ip{x,Z) denotes the deviation from the centrality. When we do iterations to minimize 
(^(x, Z) with a strict feasible start point, the first term insures (xfc, Zj^) approach the optimal 
point and the second guarantees all are feasible. More details, see [11]. A general 

outline of a potential reduction method is as follows 

Algorithm 2.1 PotentialReduction 

Input: Strictly feasible x, Z , and a tolerance e. 

Output: Updated strictly feasible x and Z . 

Repeat: 

1. Find suitable directions Sx and SZ. 

2. Find p, g G M that minimize (p{x + p5x, Z + qSZ). 

3. Update: x := x + p6x and Z := Z + q6Z. 
until duality gap < e. 

In Step 1, an obvious way to compute search directions 5x and 6Z is to apply Newton's 
method to ip. However the potential function ip is not a convex function, since the first 
term: (n + v^/r^)\og{TrFQZ + cFx) is concave in x and Z and hence contributes a negative 
semidefinite term to Hessian of ip. Wc adapt potential reduction method 2 mentioned in [11], 
which is based on the primal system only. Namely, we choose direction 5x that minimize 
a quadratic approximation of ip{x + v, Z) over all v G M™". In order to apply the Newton's 
method, the second derivative of the concave term is ignored. It is equivalent to solve the 
following linear equations: 

m 

F6ZPF + y 5x^Fi = -pFZF + F 

tl (4) 
TrFj6ZP = 0, j = l,---,m. 

we choose 6ZP as the dual search direction, sec [11]. Practically, this method seems perform 
better than method 1 which treats the primal and dual semidefinite program symmetrically. 

In Step 2, wc use plane search to choose the lengths of the steps made in the directions 
6x and 6Z, see [11]. 

n n 

minimize (n + ■uv^)log(l + cip + C2q) - ^ log(l + ppi) — ^ log(l + qvi) 

i=i i=i (^) 

S.t. Pmin ^ P ^ Pmax^ Qmin ^ Q ^ Qmax 
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Fig. 1. The objective function for plane search 



where 

~ TrF(x)Z' ~ TrF(x)Z' 

fii, . . . , and i/i, . . . , z/„ are the eigenvalues of F^^/'^6FF~^/'^ and Z^^^'^dZZ^^^'^, respec- 
tively. 

Since the objective is a quasiconvex function, we combine the Newton's direction and the 
steepest decent direction together to minimize it. If the Hessian matrix at an iteration point 
is positive definite, we use the former, otherwise we use the latter. After we get the decent 
direction 5p and 6q, we use bisection method to compute the one dimensional search of r, 
and update pk+i = Pk + r ■ 6p, qt+i = Qk + r ■ 6q, accordingly. 

Remark 2.2 Since the objective function in (5) is very sensitive at the optimizer(see Figure 
1 ), the operation must he performed carefully. A tiny numerical error can cause endless loops. 
In order to get higher accuracy result by doing more iterations, we need to set larger digits. 

Remark 2.3 We should also pay attention to choose the initial point {p,q). Usually, we 
choose (0,0) and it works well. But there comes troubles for Rump's problemfsee, 4-)- When 
n = 8, the iteration for one dimensional search does not converge. Because our methods only 
converge locally, if (0, 0) is too far from the optimizer, the iterations may not converge. If 
we denote the optimal solution as {p*,q*), according to our experiments, we find p* is always 
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attained near the largest feasible value and q* near (see Figure 1). It is always reasonable 
since when computing search direction Sx and SZ, we choose method 2, which is based on 
the primal system only. So we only compute the approximate Newton's direction of x and 
choose accessary result 5Z'p of (4) as the dual search direction. So the step length q for 5Z 
should be near 0. So we choose {p/2,0) as the initial point in plane search where p denotes 
the largest feasible value of p, and it solves the trouble of Rump's problem when n = 8. 

The algorithm PotentialReduction needs initial strictly feasible primal and dual points. 
We adapt the Big-M method when neither a strictly feasible x nor a strictly feasible Z is 
known. After modification, the primal problem becomes 



If Ml and M2 are big enough, this entails no loss of gcncrality(assuming the primal and dual 
problem are bounded). After modification, problem (6) and (7) can be written as (1) and 
(2). Then it is easy to compute its strictly feasible points, see [11] 

A brief description of the main functions contained in SDPTools is following: 

• Solution_of_Directions solves the directions 5x and 6Z from (4). 

• Expression_of_PlaneSearch gets the objective of (5). 

• StepsLength_PlaneSearch solves problem (5) and get step length p,q for x,Z respec- 



• BigM_Case3 transforms (6), (7) into (1), (2) respectively, and computes a strict feasible 

start point {xq,Zq). 

• Solve_SDP_Method2 solves the SDP (1) and (2). 

3. Certified Global Optimum of Rational Functions 

In SDPTools, we apply the SDP algorithm (2.1) to compute and certify the lower bounds 
of rational functions with rational coefficients. 



minimize (Fx + Mi t 

s.t F{x) + tl>0, 
TrF < M2, 
t > 



(6) 



The dual becomes 



maximize —TrFo{Z — zl) — M2Z 
s.t TrFi{Z - zl) = Ci, i 
TrZ < Ml , 

Z > 0,z> 



m. 



(7) 



tively. 




(where g{^) > for ah ^ G M") 



(8) 
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where f{C),g{6.) G Q"- The number r is a lower bound of (8) if and only if the polynomial 
f{x) — rg{x) is nonnegative. Therefore we focus on the following minimization of a rational 
function by SOS: 



: sup r 
reR,W 



S.t f{x) - rg{x) = md{x)'^ ■ W ■ md{x) 



(9) 



where md{x) is the column vector of all terms in Xi, . . . , X„ up to degree d. A detailed de- 
scription of SOS relaxations and its dual problems, sec [5]. The problem (9) is a semidefinite 
program and can be written as (2) . With packages for solving SDP, we can obtain a numerical 
positive semidefinite matrix W and floating point number r* which satisfy approximately: 

f{x) - r*g{x) « md{xY ■ W ■ md{x), W^O (10) 

To certify r*, we convert r*,W to rational ones and project orthogonally W onto the affine 
linear hyperplane: 

X = {^1^"^ = A, f{x) — rg{x) = md{x)'^ ■ A ■ rridix), for some r} (11) 

and hope W to be positive semidefinte after projection. 

However, if we run the fixed precision SDP solver in Matlab, the output W are too coarse 
to be projected into the cone of positive semidefinite matrices. In [4] they refined the r*,W 
by Gauss-Newton method. Here our high precision SDP solver shows its advantage. It is 
implemented in Maple, which has arbitrarily high precision, so it can compute r*,W with 
high accuracy. Without refinement, the projection can be done successfully. 

Usually, it is not easy to find a strictly feasible point for (9) and we need the Big-M 
method to avoid it. After convert (9) to the form (7), the SOS relaxation of (8) becomes 



r := sup r — M2Z 

reM.,W 

S.t. f{x) - rg{x) + z{md{xY ■ md{x)) = md{xf ■ W ■ md{x), 
wyo, = W, z>0 



(12) 



It is obvious to see that problem (9) and (12) are equivalent. And the variable z is required 
to be at the optimizer of problem(12). We define 

X = {A \ = A, f{x) - rg{x) + z{md{x)'^ ■ md{x)) ^^^^ 
= m,d{x)'^ ■ A ■ md{x), for some r, z} 

Note that x will meet x a-t the optimizer (W*,r*), see Figure 2. 

After each iteration, wc get a point (r, W) on x- Then we convert (r, W) to rational 
ones and project matrix W onto x, and denote the nearest matrix in x by W. Because 
the problem(9) and (12) are equivalent, we can expect that after the first few ones, at each 
iteration we can obtain a positive semidefinite matrix W and a certified r. Different from 
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Fig. 2. Rationalization of SOS 



the method in [4] which only certifies a given r, we can get series of certified lower bounds 
rn- The more iterations we do, the better r„ we get. The process is shown in Figure 2. 
In order to reduce the problem size, we also exploit sparsity. Given a polynomial p{x) = 
the cage of p, is the convex hull of sup(p) = {a\ pa 7^ 0}. 

Theorem 3.1 ([6]) For a polynomial p, C{p'^) = 2C{p); for any positive semidefinite poly- 
nomial f and g, C{f) C C{f + g); if f = q] then C{gj) C lc(/). 

With the above theorem, we can remove the redundant monomials in (9). 

SDPTools also has functions to compute the convex hull of given finite points in n dimen- 
sional vector space. They are mainly based on the Quickhull algorithm in [1] which assumes 
the affine dimension of the input points is also n. If the affine dimension, denoted by r, is 
less than n, SDPTools provides a function to compute r and affinely transforms the input 
points to a r dimensional vector space where we can use the Quickhull algorithm. For a 
given point, SDPTools has a function to judge whether it is in the convex hull. 

SDPTools provides a function to compute the monomials appearing in (9). For example, 
when solving Rump's problem, the monomials appearing in mrf(x) are obtained by proof in 
[4] while we can find them automatically by SDPTools! 

The followings are the main functions to compute and certify the lower bound of a given 
rational function: 

• affineDims computes the affine dimension of input points. 

• affineTrans affinely transforms the given points to a r dimensional vector space. 

• convexHull computes the convex hull of given points. 

• inConvexHull judges whether p is in the convex hull defined by the points in 5. 

• getSDP transforms problem (9) as form (2). 

• projSOS converts r, W to rational ones and projects W onto x- 

• certifiSOS computes and certifies the lower bound of an input rational function. 
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4 X 15 


1.76 


0.00028973187527968191 


0.00028973187527968193 


7 


1 


75 


5 X 15 


11.36 


0.00003418506980008284 


0.00003418506980008285 
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2 


75 


5 X 15 


12.49 


0.00000390543564975572 


0.00000390543564975573 


9 


1 


75 


5 X 15 


84.12 


0.43600165391810484612e-06 


0.43600165391810484613e-06 


10 
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75 


5 X 15 


92.79 


0. 478393956877097593260-07 


0. 478393956877097593270-07 
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5 X 15 


622. 0:i 


0.51787 19097 1 l()99():)33()c-()8 


0.51787 19097 1 l()99()83:!lo-()8 
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5 X 15 
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0.55458818311631347611e-09 


0.55458818311631347612C-09 
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5 X 15 


3800.0 


0.58866880811866093129e-10 


0.58866880811866093130e-10 


14 
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100 


5 X 15 


3800.0 


0.62024449920539050219e-ll 


0.62024449920539050220e-ll 


15 


1 


120 


6 X 15 


15000 


G.64943654185809512879e-12 


0.64943654185809512880e-12 


16 


2 


120 


6 X 15 


23000 


0.67636042558221379057e-13 


0.67636042558221379058e-13 



Table 1. The certified lower bounds 



4. A numerical experiment: solving Rump's Model Problem 

The following introduction to Rump's model problem (see, [7]) is mostly from [4]. This 

problem is related to structured condition numbers of Toeplitz matrices and polynomial 
factor coefficient bounds, asks for n = 1, 2, 3, ... to compute the global minima 



\\PQ\\ 



2 



n n 

s. t. p{z) = Y,Pi^'~\Q{^) = ^ \ io}. 

i=l i=l 

It has been shown in [8] that polynomials P, Q realizing the polynomials achieving jji^ must 
be symmetric (self-reciprocal) or skew-symmetric. Thus the problem can be rewritten into 
three optimization problems with three different constraints 



k = 1 

k = 2 
k = 3 



Pn+i-i = Pi, qn+i-i = qi, 1 < « < n, 

Pn+l-i = Pi, Qn+l-i = -Qi, I <i <n, 
Pn+l-i = -Pi, Qn+l-i = -Qi, I <i <n, 



and the smallest of three minima is equal to /in- For all three cases, we minimize the rational 
function f{X)/g{X) with 

2n n n 

f{X) = \\PQ\\l = J2{J2 Pi^j)'^ 9{X) = \\P\\l\\Q\\l = {J2phiJ2Q]) 

k=2 i-\-j=k i=l j=l 

and the variables X = {pi, . . . ,Pn(p)} U {qi, . . . ,(?„(q)}, where n{P) = n{Q) = \n/2\. 

In this paper, we use higher precision SDP solver in SDPTools to solve Rump's model 
problem and obtain much better certified lower bounds, see Table 1. 



5. Conclusion 

SDPTools is a Maple package having the following functions: to solve a general SDP, to 
compute and certify the lower bounds of rational functions. It is also a tool to compute the 
convex hull of given finite points in order to explore the sparsity structure. 
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SDPTools is still under development, we hope to implement more efHcient algorithms to 
solve SDP, to give upper bounds for Mi, M2 when applying Big-M method in SOS relaxation 
(12), to detect the infeasibility for SDP and to explore more sparsity structures for problems 
having large size. 
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